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Abstract 

We propose a means for constructing highly accurate equations 
of state (EOS) for elemental solids and liquids essentially from first 
principles, based upon a particular decomposition of the underlying 
condensed matter Hamiltonian for the nuclei and electrons. We also 
point out that at low pressures the neglect of anharmonic and electron- 
phonon terms, both contained in this formalism, results in errors of 
less than 5% in the thermal parts of the thermodynamic functions. 
Then we explicitly display the forms of the remaining terms in the 
EOS, commenting on the use of experiment and electronic structure 
theory to evaluate them. We also construct an EOS for Aluminum 
and compare the resulting Hugoniot with data up to 5 Mbar, both 
to illustrate our method and to see whether the approximation of 
neglecting anharmonicity et al. remains viable to such high pressures. 
We find a level of agreement with experiment that is consistent with 
the low-pressure results. 

1 Introduction 

Over approximately the last sixty years, numerous models and techniques 
have been developed for creating equations of state (EOS) for a variety of 
materials that are valid up to very extreme pressures (tens of Mbar) and 
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temperatures (several eV). In the EOS community at the national labora- 
tories, for instance, we have often used models based on the Mie-Griineisen 
EOS together with the Thomas- Fermi or Thomas- Fermi-Dirac model (or one 
of its modifications) to include the contributions from the electrons (see 
for examples). The models usually contain enough independent parameters 
to adjust the EOS until it correctly reproduces the experimentally measured 
Hugoniot (and perhaps a few other data points), but it is generally an open 
question how accurate the EOS is away from the Hugoniot. In this paper we 
argue that for one class of materials, elemental solids and liquids, our under- 
standing of the underlying condensed matter Hamiltonian for the nuclei and 
electrons has grown to the point that we can construct highly accurate EOS 
from essentially first principles, and we also propose a means for doing so. 
We also argue that, since the underlying physics is well understood, an EOS 
derived this way should have the right functional form, even if we are unsure 
of the values of some of its parameters; thus, if the resulting EOS is shown 
to be accurate in one thermodynamic region (say, along the Hugoniot), then 
we can be confident that it is roughly equally accurate elsewhere. 

In this formalism, the EOS in the solid phase depends on a decomposition 
of the Hamiltonian due to Wallace (see Chapter 1 of ) , extending the work 
of Born [3] to metals as well as insulators; the resulting free energy contains 
terms describing the harmonic motion of the nuclei about their lattice sites 
(phonons), thermal excitation of the electrons from their ground state, an- 
harmonic corrections to the nuclear motion (represented as phonon-phonon 
interactions), and interactions between the electron excitations and the nu- 
clear motion, represented as electron-phonon interactions. (Please note that 
this description is exact; all of the physics contained in the true Hamilto- 
nian of the system is included here. Specific EOS models usually neglect the 
anharmonic and electron-phonon terms, arguing that anharmonicity is small 
and making reference to some form of the Born-Oppenheimer approximation; 
we will take a somewhat different route, commenting on approximations be- 
low.) A recently developed theory of the dynamics of monatomic liquids (see 
[3] for a review) uses the same Hamiltonian to derive a liquid free energy 
which is quite similar to the expression for a solid, with additional terms 
accounting for the fact that the liquid, as opposed to the solid, traverses 
many potential valleys and thus sees the boundaries between them. For both 
phases, the resulting free energies have been compared with experimental 
data in the low-pressure regime (P < 100 kbar), with the following results 
(Sections 17-19 and 23 of 0): 
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(a) Molecular dynamics (MD) calculations of the anharmonic contribu- 
tion to the entropy of several solids match experimental entropy data to the 
accuracy of the data themselves. 

(b) Low-temperature (T < 20 K) calculations of the electron-phonon term 
for several solids lead to predictions that also match experimental entropy to 
the accuracy of the data. 

(c) Theoretical arguments show that the electron-phonon contribution is 
entirely negligible except when the electronic contribution dominates the free 
energy, such as in metallic solids at low temperatures. 

(d) For the 27 elemental solids for which accurate data are available from 
low T (but not too low; see point (c)) to the melting temperature T m , the 
free energy excluding the anharmonic and electron-phonon terms accounts 
for the experimental thermal energy and entropy to an accuracy of 5% (in 
fact, an accuracy of 2% for all but about five materials). 

(e) For the 6 elements in the liquid phase for which accurate data are 
available at temperatures up to around 3T m , the effect of neglecting the 
anharmonic, boundary, and electron-phonon contributions to the energy and 
entropy is similarly small. 

This tells us that at low pressures, we can neglect the anharmonic, bound- 
ary, and electron-phonon terms in both the solid and liquid free energy (which 
happen to be the hardest terms to calculate), and the resulting thermal en- 
ergy and entropy are both simple in form and accurate at the 5% level. It 
is for this reason, not an appeal to the Born-Oppenheimer or other approxi- 
mations, that we know we can simplify our EOS and what the results of the 
simplification will be, at least at low pressures. 

In this paper we do two things: (1) We describe in more detail this 
framework for constructing EOS and discuss the theoretical and experimental 
inputs needed to implement it, and (2) we construct a sample EOS, neglecting 
anharmonic, boundary, and electron-phonon terms, both to illustrate the 
method and to discover whether points (d) and (e) above continue to hold 
in the high-pressure regime. We use Aluminum as our sample because of the 
availability of extensive electronic structure calculations, up to a compression 
of three, and highly accurate shock Hugoniot data, which provide a test of our 
EOS through both phases to pressures of around 5 Mbar. In Subsection 12. II 
we develop the general theory of the solid EOS, and in Subsection 12 . 21 we do 
the same for the liquid. In Subsection 13.11 we construct our sample EOS for 
Al, comparing it with other EOS work, and in Subsection 13 . 21 we compute the 
Hugoniot predicted by the EOS and compare it with experimental data. The 
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results are encouraging. Finally, we review our work, discuss the advantages 
and disadvantages of this formalism (and how to address the disadvantages), 
and suggest directions for future development. 

2 General theory 
2.1 Solid phase 

The condensed matter Hamiltonian, decomposed as described above, consists 
of terms describing the motion of the nuclei in a potential generated by the 
electrons in their ground state, plus additional terms that lead to the thermal 
excitation of the electrons and describe their interactions with the nuclear 
motion. With this Hamiltonian, the Helmholtz free energy per atom for a 
solid at temperature T with volume V per atom takes the form 



F*(V, T) = $ S (V) + F*(V, T) + F^V, T) + F^V, T) + F^V, T). (1) 



Here <3>q is the static lattice potential (the electronic ground state energy 
when the nuclei are fixed at their lattice sites); it depends on the particular 
crystal structure. F^ h is the contribution from the harmonic motion of the 
nuclei about their lattice sites, F* x represents the thermal excitation of the 
electrons when the nuclei are fixed at their lattice sites, F^ nh accounts for the 
anharmonicity of the nuclear motion (which may be represented as phonon- 
phonon interactions), and F e s p expresses the interactions between the electron 
excitations and the nuclear motion, represented as electron-phonon interac- 
tions. (We emphasize again that this free energy is exact; it includes all of 
the physics present in the Hamiltonian.) The discussion in the Introduction 
justifies our approximating the solid free energy as 



so let us now consider the forms of F* h and and the parameters on which 
they depend. The phonon term in the Hamiltonian describes harmonic mo- 
tion, which leads uniquely to the free energy of lattice dynamics: 



where j3 = 1/kT and g s (uj) is the distribution of phonon frequencies in the 
Brillouin zone. (Note that g s (u) is volume dependent.) Sometimes we require 



F S (V, T) = $ S (V) + F°(V, T) + F^V, T), 



(2) 




(3) 
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not the full Eq. (J3J) but only its high- and low-temperature limits, for which 
we need not the full g s (uj) but only three of its moments, expressed in terms 
of the characteristic temperatures 0q, 6^, and 9| defined by 



lnke s 
kQ{ 
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o(Mbz 



tt((M 2 )bz 



1/2 
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where (• ■ -)bz indicates an average over all the frequencies in the Brillouin 
zone. Then the following limits hold: 
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as 



T —> 



and 



F*(V,T) = -3kT 



1 f T \ 1 f Q 2 



+ 



at high T. 



(5) 



(6) 



The leading term in Eq. (jHJ) describes purely classical nuclear motion, while 
the series of terms in powers of T~ 2 are quantum corrections. Keeping only 
the first quantum correction, the thermodynamic functions derived from Eq. 
© are accurate to 1% at temperatures above 

The electronic excitation free energy can be expressed generally as an 
integral function of the electronic density of states per atom, n s (e), and the 
Fermi distribution ^ 

where (3 is still 1 / kT and \x is the chemical potential. If each atom contributes 
Z electrons to the valence bands (notice that Z is not necessarily the atomic 
number), with the lowest valence energy set to zero, then \i is a function of 
T determined by the normalization condition 



/ n\e)f(e)de = Z. 
Jo 

The electronic free energy is then 



(8) 



Ki(y,t) 

jiZ — 



en s (e) de - kT n s (e) ln[l + e^"^] dc, 
Jo 



(9) 
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where ep, the Fermi energy, is the value of // when T = 0. The second term 
on the right hand side of Eq. Q is the subtraction of the electronic ground 
state energy, which ensures that — > as T — ► 0. This property makes 
sense if represents purely thermal excitation of the electrons. (It also 
avoids double counting of the energy, as the electronic ground state energy 
is already represented as $q.) 

We see from this discussion that to evaluate the terms in Eq. (J2J) for 
the solid free energy we require three unknown functions: $q, 5 |S ( a; ) ( or ®o> 
0f, and 0| if we are concerned only with the high- and low-T limits), and 
n s (e) (and the associated quantities Z and ep). These can be determined in 
various ways: compressibility data and diamond anvil cell data can be used 
to construct $q(V); neutron scattering experiments can determine g s (w) or 
its various moments at P = 1 bar; and for many elements all three of these 
functions can be computed reliably using electronic structure theory. (Or 
one could use results from multiple sources in combination, which is often 
an option with $q and is basically a necessity with g s (u>).) One must keep 
in mind, however, that the accuracy of one's answers will be limited by the 
accuracy and range of applicability of these functions, regardless of how they 
are determined. 

2.2 Liquid phase and two-phase region 

According to the theory of liquid dynamics reviewed in [3] , the same Hamil- 
tonian that gave us the solid free energy leads to a similar form for the free 
energy of a monatomic liquid. In this theory, the region of the many-body 
potential surface in which the system moves in the liquid phase is dominated 
by a large number of intersecting nearly-harmonic valleys, called "random" 
valleys because they correspond to particle configurations which retain no 
remnant crystal symmetry, and which are all macroscopically identical. In 
particular, the valleys all have the same distribution of normal mode fre- 
quencies, and they all have the same depth (which, as in the solid case, is 
the electronic ground state energy when the nuclei are fixed at the valley 
minimum). The resulting liquid free energy per atom is 



All of the terms correspond to their solid counterparts with the following 
exceptions: 



F\V,T) 



TO + *Sh(V, T ) + Ki(V, T) + Fl h (V, T) + 
FUV,T)-kT\nw. 



(10) 
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(1) $q, now called the static structure potential, is the depth of a typical 
valley in which the liquid system moves. 

(2) The normal mode spectrum appearing in is that of a typical liquid 
potential valley, not the unique solid potential valley. 

(3) The term F^ h includes corrections due both to anharmonicity and 
to the fact that the potential valleys have boundaries, which the liquid (as 
opposed to the solid) encounters as it transits from valley to valley. 

(4) The extra term — kT In w corresponds to an increase in entropy of 
fclnw per atom; the value lnw » 0.8 is estimated from entropy of melting 
data of the elements (again, see [I] for details). In liquid dynamics theory, 
this term is due to the hypothesis that the number of potential valleys among 
which the liquid moves is of order w , where N is the number of atoms in 
the system. 

We emphasize that the same Hamiltonian gives rise to both Eq. (JTJ and 
(jTUJl : the differences are that the potential is expanded about different equi- 
librium configurations in the two cases, and that the region of configuration 
space over which the liquid moves is obviously far larger than the space 
available to the solid (hence the — kT In w term). 

Again making the approximations discussed in the Introduction, our form 
for the liquid free energy becomes 

F\V,T) = & (V) + Fl h (V,T) + F* l (V,T) - kTlnw, (11) 

and the additional term — kTlnw is fully determined by setting Inw = 0.8, 
as mentioned above. The form of the phonon term is dictated by a central 
hypothesis of liquid dynamics theory: The motion of the liquid consists of 
oscillations in the macroscopically similar valleys described above together 
with occasional transits between valleys; the transits are of such short du- 
ration that they do not contribute to the thermodynamics to lowest order. 
Thus we will take FL to have the same form as the solid phonon term, Eq. 
Q, with a possibly different phonon frequency distribution g {oS). The elec- 
tronic excitation term for the solid was derived using only the assumption 
that the electrons are thermally distributed over the available states using 
Fermi statistics; all of the information about the configuration of the nuclei is 
contained in the density of states. Hence F^ also takes the same form as the 
corresponding solid term, Eq. Q, with a density of states n l (e) appropriate 
for the liquid phase. (What this means is discussed briefly below.) 

The liquid and solid EOS together determine the melting temperature as 
a function of pressure T m (P) by the requirement that the solid and liquid 



7 



Gibbs free energies match along the melt curve, or 



G s (P,T m (P)) = G\P,T m (P))- (12) 

Once the melt curve is determined, one can use the solid and liquid EOS 
separately to compute V^(T) and V^(T), the atomic volumes of the solid 
and liquid at melt as functions of temperature. (Of course, using T m (P) one 
can express and as functions of pressure equally well.) 

In the case V^(T) < V^JT) for all T, which we assume here, the compu- 
tation of the full two-phase EOS (F, E, S, and P) for any V and T proceeds 
as follows. If V < V^(T) for the given T, the system is in the solid phase; F s 
is computed as described in Subsection 12.11 and the other functions follow 
from 

' - -(£),■ 

- - -(£),■ 

If V > V^(T), the system is in the liquid phase; F 1 is computed as described 
in this Subsection, and the other functions follow from expressions analogous 
to Eqs. (JT3J). If V^(T) < V < V^(T), the system is in the two-phase region; 
defining rj, the volume fraction of the system in the liquid phase, by 

V - V s (T) 

V Vi(T)-V^(TY { ] 

the thermodynamic functions are 

F(V,T) = r / F z (^(T),T) + (l-r ? )F s (^(T),T), 

E(V,T) = V E l (Vi(T),T) + (l- V )E s (V^(T),T), 

S(V } T) = v S l (Vi(T),T) + (l- v )S s (V^(T),T) } 
P(VT , F l (Vj(T),T)-F%V^T),T) 

m ] ~ vm-vi(T) • (15) 

Finally, we note that just as with the solid, to evaluate the terms in Eq. 
(|TT|) we require three unknown functions: $q, g l {oS) (or Oq and ©2, since we're 
not likely to need the low-T limit), and n l (e). In contrast to the solid case, 
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these functions are generally not available experimentally. (It is possible that 
one might be able to compute $q using liquid compressibility data, but we 
suspect that this will be very difficult.) However, for many materials these 
functions should be computable using electronic structure theory, proceeding 
much as one would in the solid case except that the nuclei would be arranged 
not in a crystal configuration but in a disordered structure characteristic of 
a "random" valley in the liquid potential surface jlj. To our knowledge very 
few such calculations have been attempted; the only ones we are aware of 
are <3>q and g l {ui) at a single volume for liquid sodium in [5] (the results are 
referred to in jl] and a graph of g l {oj) using their results appears as Fig. 1 in 
[U]). Another function that is sometimes available is the melt curve T m (P), 
but this curve cannot be chosen independently of the others, since the solid 
and liquid EOS determine it jointly; but this can be an advantage, since if 
T m (P) is known from experiment, for example, it can be used to compute 
one of the other needed functions if it is not otherwise available. In fact, this 
is how we will determine $g in our example EOS, to which we now turn. 

3 An example: Aluminum 

To illustrate the application of the theory we've described, we will now con- 
struct an EOS for Aluminum, which has been the subject of extensive elec- 
tronic structure calculations and for which a great deal of high-quality ex- 
perimental data are available. We will then compare the Hugoniot predicted 
by our EOS with data up to pressures of approximately 5 Mbar; this will tell 
us whether the approximations we discussed in the Introduction (neglecting 
anharmonic, boundary, and electron-phonon effects), known to be accurate 
at low pressures, continue to be reasonable in the high-pressure domain. 

3.1 Constructing the EOS 

We recall from Subsection 12.11 that the solid EOS requires three functions: 
$q, g s (u>), and n s (e). Since we will be testing the EOS by comparison with 
Hugoniot data, we will always be in the high-T region (except for one brief 
low-T excursion; see below), so we use Eq. © for F^ h instead of Eq. (jSJ); this 
means that we require only Oq, O^, and 0| in place of g s (ui). 

To determine these functions, we began by consulting the results of den- 
sity functional theory (DFT) calculations carried out in the local density ap- 
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proximation by Straub et al. jjj. They worked with fee and bee Al at atomic 
volumes from 37 a 3 , to 160 ajj, where do is the Bohr radius, corresponding 
to densities from 8.17 g/cm 3 to 1.89 g/cm 3 (the density of Al at 293 K and 
1 bar is 2.700 g/cm 3 ). Their calculations indicate a T = transition from 
fee to bee at 51 a 3 ,, corresponding to p = 5.93 g/cm 3 , but we will neglect 
this phase change and treat solid Al as an fee crystal for two reasons: The 
DFT calculations themselves suggest that the effect of the phase change on 
the thermodynamic functions will be quite small; and we know from experi- 
ment that the solid-liquid transition on the Hugoniot takes place well before 
reaching the density of concern (see Subsection 13.2)1 . so we are confident of 
our assumption of fee along the Hugoniot until melting. However, this as- 
sumption may have an effect on the liquid EOS at high densities, which we 
will comment on below. (Other electronic structure work, discussed on pp. 
89-90 of |H] , suggests the possibility of an hep phase between the fee and bee 
phases, but as jH] also mentions, no experimental signature of this phase has 
been found, so we will proceed under the assumption of a single solid phase.) 

Straub et al. computed $g for fee by fitting their results to a Birch- 
Murnaghan form, 

with coefficients 

V b = 106.302 a 3 , c = -287.7832 mRy, 

c 2 = 761.2029 GPa, c 3 = 1319.036 GPa, 

c 4 = -13 661.06 GPa, c 5 = 50 315.53 GPa. (17) 

The were determined by computing the bulk modulus and four zone- 
boundary phonons at several volumes, and these results were used to calibrate 
a pseudopotential model at each volume. The pseudopotential was then 
used to calculate phonon frequencies throughout the Brillouin zone, allowing 
the determination of ©q, Q\, and 0|. Their results are shown in Tabled 
and Figure [T] (The full set of results was not reported in (7j.) To check 
these results, Straub et al. compared experimental phonon moments for Al at 
T = 80 K and P = 1 bar based on Born- von Karmen fits to neutron scattering 
data 9j with their predictions interpolated to the appropriate atomic volume 
of 110.7 Oq- The experimental points, also shown in Figure ^ are in very 
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V (a§) 


©o (K) 


©i (K) 


©1 (K) 


111.97 


278.09 


386.55 


387.20 


106.65 


304.63 


423.81 


424.86 


93.318 


381.43 


532.00 


534.48 


74.655 


525.01 


735.49 


741.68 


55.991 


741.62 


1044.7 


1058.3 


37.327 


1109.5 


1575.0 


1605.4 



Table 1: DFT calculations of G s , G\, and G s 2 from [7j. 




Figure 1: Qq, Q^, and Q s 2 as functions of atomic volume from the DFT 
calculations in j7] and experimental data from jHj. Our functional forms are 
also shown. 
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good agreement with their calculations; hence these results for the 9^ are 
acceptable for use in our EOS without modification. To determine the 0^ at 
any volume, we first constructed a functional fit to the Oq points, with the 
result 

17319 9 

Q a Q (V) = 2852.69 + — + 2.33667 V - 633.858 ln(V), (18) 

where Q s is in K and V is in a^. Then we noted that according to the DFT 
results both Q\ and B| approximately equal e 1//3 0q, so we made the same 
approximation using Eq. ()18|) for 0q to determine and 0| at any volume. 
These functions are also shown in Figure ^ 

The DFT calculations also provided data on the electronic density of 
states n s (e). Graphs of n s (e) for fee and bec Al at atomic volume 112.0 cIq 
(corresponding to P = and T = 295 K) are shown in Figure along with 
the free electron n s (e), for which 

at V = 112.0 Oq and Z = 3. The Figure shows that the free electron model is 
a good approximation for either crystal structure, for electronic excitations 
to around | Ry. The same match, at a slightly poorer level of approximation 
and for excitations to around 1 Ry, is found at our smallest atomic volume 
of 37 Oq. For all volumes of our study and temperatures up to T m , the 
total electronic excitation contribution to the energy, entropy, and pressure 
is at most 5%, so the error introduced by using the free electron n s (e) in 
our calculations is negligible. Making this approximation, the normalization 
condition from Eq. (jHJ) becomes 



2 



*i/a(&0 = ^fY /2 , (20) 

which determines function of V and T, and Eq. Q for the free energy 

becomes 

w^4-^ r 3i-H (2i) 

where Z = 3. The F n (x) are the standard Fermi integrals; their properties are 
discussed on pp. 332-334 of [TU] and their values for various n are tabulated 
in [HI EES- 
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Figure 2: n s (e) for bcc and fee Al at atomic volume 112.0 from the cal- 
culations in The free electron n s (e) at this volume and Z = 3 from Eq. 
fll9j) is shown for comparison. (From 



The solid EOS that results from assembling all of these functions is re- 
liable over a large range of volumes and temperatures; however, it is not in 
perfect agreement with the highly accurate experimental data that are avail- 
able at low pressures. Specifically, experiments on Al at T = and P = 
show that [7] 

V = 110.6 al, E = -249 mRy, 
d R 

£? = 79.4GPa, ^ = 4.7, (22) 

but the EOS yields Vo = 107.3 Oq, which is outside the experimental error. 
Therefore, we chose to make a small correction to our purely theoretical free 
energy to agree with experiment. These quantities are obviously determined 
by the T = form of the free energy, Fq = $ + ^kQ\ (see Eq. (jSJ); since 
is already in excellent agreement with available experiment, we chose to 
modify only $q. To proceed, we noted that the data determine P , the T = 
pressure, in the vicinity of V = Vo by the relation 

(V - V ) 2 + • • • 

Vo 



Po(V) = Po(V ) + 



dPo 

dV 



(V - Vo) + - 



1 d 2 P n 



V'n 



2 dV 2 
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-t^-^ + w( 1+ ^ 1(K -^ + - (23) 



while at higher compressions we have no information to supplement the elec- 
tronic structure result; so we decided to construct a $q that correctly re- 
produces Eq. (J23J) near Vq but smoothly interpolates to Eq. (fTH|) at lower 
volumes. To do this, we computed Pq at 10 volumes between 110 and 
111.25 do using Eq. (J23J), and we also computed P = —dF^/dV using the 
above form for Fq, with $q from Eq. (|16p. at 23 volumes between 30 a,Q and 
41 Oq. We then performed a least-squares fit to these points using an expres- 
sion similar to the Birch-Murnaghan form, but carried to a slightly higher 
order; after integrating, adjusting the constant of integration to correctly 
match E from Eq. and subtracting off f&O®, we had a new $q given 
by 

_ /TrN , 2.07608 x 10 7 4.61515 x 10 s 

$ S (V) = -1.64615 x 10 6 + - 



^2/3 y4/3 

5.71249 x 10 9 5.49998 x 10 10 3.71978 x 10 11 

1 V* ^873 1 yw/3 

1.66284 x 10 12 4.41118 x 10 12 5.25064 x 10 12 
+ 



V 4 V 14 / 3 V 16 / 3 

+ 2.26789 x 10 7 _ ^ ^ y + ^ ^ ^ ^ 

where $q i s i n m Ry and is in a^. This which reproduces the data in 
Eq. (}2*2*|) and interpolates smoothly to the DFT curve at higher compressions, 
is what we use in our EOS instead of Eq. ([16)1 . The T = pressure- volume 
curves constructed using both the original and new $q are shown in Figure 

El 

Our choice of a Birch-Murnaghan-like form was dictated by the fact that 
the Straub et al. result provides most of our information about the shape 
of $q; so our goal was to preserve that form insofar as was possible, in- 
terpolating back to their result as quickly as we could without introducing 
enough curvature to compromise agreement with dB /dP. This correction 
to $o constitutes a small change to the overall EOS; the effect of this change 
on the Hugoniot will be considered in the next Section. This modification 
completes the full solid free energy, so we can now consider the liquid. 

For the liquid we need the same three functions that we needed for the 
solid, and we must also consider the melting curve T m (P). Having chosen to 
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Figure 3: The T = pressure-volume relations calculated using the original 
$q and the new $q we constructed. Notice how they differ in the vicinity of 
V = 110.6 al but then agree more closely at lower volumes. 
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use Eq. (jUJ) for F* h , we certainly did the same for F^ h , since the Hugoniot 
will obviously enter the liquid region only at rather high temperatures; thus 
we needed only O and 2 . From experiment we know that Al is what is 
called in liquid dynamics a "normal melting element" (the entropy of melting 
at constant density is approximately 0.8 k per atom), and we argue in [3] 
that O in solid and liquid phases of such an element are approximately 
equal. (Experimental and computational work supporting this conjecture 
are also discussed in 0j.) Thus we took the liquid to have the same 0o as 
the solid. It is also true that in the liquid, T is typically much larger than 
G| (for example, in liquid Al at normal density T > 20|), rendering the 
first quantum correction to Fp h negligible (roughly 1% at normal density), so 
even if 2 differs from 02 by 25% or more, the impact on the phonon term 
will be very small; therefore we also used the same 02 in the liquid as in the 
solid. 

Since the free electron model approximates the DFT result for n s (e) so 
well for both fee and bec Al (Figure|2J), which correspond to two valleys in the 
many-body potential surface with rather different structures, we also expect 
this model to be a good approximation for n l (e), the density of states for the 
structure characteristic of a liquid. Since at all volumes and temperatures 
up to 5T m (the relevance of this number will appear below), the electronic 
contribution to the thermodynamic functions does not exceed 25%, the error 
introduced by the free electron model is still acceptable. 

We fixed the only remaining term in Eq. (|11J1. $q, by the requirement that 
the Gibbs free energies of the solid and liquid match along the Al melting 
curve. Boehler and Ross [T3] suggested that 



on the basis of their experimental work up to 80 GPa (0.8 Mbar), and exper- 
iments by McQueen et al. and Hanstrom and Lazor ^2] and theoretical 
work by Pelissier [TH] suggest that this curve continues to be accurate up to 
200 GPa. In the absence of evidence to the contrary, we took Eq. (|23j) to 
be valid to higher pressures as needed. (As we will see later, our EOS will 
assume Eq. (|25|) no higher than 400 GPa.) We computed $q as follows: We 
made a guess for $(, not very different from $q, and then we used it and Eq. 
f!25|) to calculate the difference between the two Gibbs free energies, 




(25) 



AG(P) = G s (P,T m (P)) - G l (P,T m (P)), 



(26) 
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at several hundred values of P over the entire pressure range considered in 
this study. We also calculated the liquid melt volume V^(P) at each P. If the 
rms average of Eq. ()26|) over all calculated points was not sufficiently small, we 
used the following easily verified fact: To first order, a small change 5& pro- 
duces a small change 5G l (P,T m (P)) given by 5G l (P,T m (P)) = 5$ (V^(P)). 
Thus we performed the substitution 

^ l (V)^^ (V)+AG(Pi l (V)), (27) 

where AG was computed by Eq. and P}n{V) is the inverse of V^(P), and 
calculated Eq. (}2T)j) again. We iterated until the rms average was sufficiently 
small (less than 0.001 mRy in our case), giving us the needed $(,, which is 
shown in Figure 0] along with $q. We recorded as a list of points, and 
we did not create a functional fit for it; instead we used an interpolator to 
calculate it and its derivative when needed. 

It is at this point that the existence of other solid phases in Al, discussed 
earlier, affects the EOS of the liquid. It is likely that the liquid region borders 
the fee crystal only over part of its boundary, beyond which the liquid borders 
the bec region or other solid phases. In such a case, at sufficiently high 
pressures we should use the free energy appropriate for that solid phase, not 
the fee free energy, in Eq. (}2l)j) . This suggests that $q may become inaccurate 
beyond densities in the neighborhood of 6 g/cm 3 , where the T = fcc-bcc 
phase transition occurs. We will take this fact into consideration when we 
discuss the limits of applicability of the EOS below. 

Once we had the full solid and liquid EOS, we then solved Eq. (fT^j) directly 
to compute T m (P), verifying that we had reproduced the Boehler-Ross curve; 
our result is shown in FigureEl together with the data from [T^ ITo] and some 
points from Pelissier's theoretical curve. (According to [T3|, their data point 
at 125 GPa marks the onset of melting along the Hugoniot. We will comment 
on this in the next Subsection.) 

Now that we have the full two-phase EOS, it is profitable to compare our 
work with another EOS for Al, due to Moriarty et al. [TZj. These authors 
also use a full lattice dynamics treatment of the crystal phonons, although 
they calculate their g s (uj) two separate ways, using both Moriarty's General- 
ized Pseudopotential Theory (GPT) and a local pseudopotential model with 
parameters chosen to match solid-phase EOS data. We strongly prefer to 
rely on DFT results, as we believe DFT has reached such a stage of maturity 
that it more accurately captures the physics contained in the real Hamilto- 
nian of the system, which as we have emphasized we believe to be understood 
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Figure 4: $q determined by matching the liquid and solid Gibbs free energies 
along the melt curve. $q is also shown for comparison. 
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Figure 5: The melt curve T m (P) computed from our full solid and liquid 
EOS (which reproduces the Boehler-Ross curve, Eq. (|25|)). the experimental 
data from ^J^|, and points from the theoretical curve in [To] . 
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in sufficient depth that it should underlie all of our work. Second, in their 
treatment of the liquid phase Moriarty et al. rely on fluid variational theory, 
described in detail in [15] . to compute the least upper bound to the "real" 
liquid free energy (from a liquid Hamiltonian based on pseudopotentials) 
that can be obtained from the free energy of a reference system; Moriarty 
et al. investigate hard-sphere, soft-sphere, and one-component plasma refer- 
ence systems before settling on the soft-sphere system as providing the best 
bound. Based on the investigations summarized in jl], we claim that we 
have the actual Hamiltonian of the liquid itself, not a Hamiltonian based on 
pseudopotential theory; furthermore, this Hamiltonian decomposes naturally 
into a dominant term that produces a free energy that can be used directly 
(instead of requiring approximation by the free energy of a reference system) 
and remaining terms whose contributions to the free energy are known to be 
small (see the Introduction). The same point we made above for the solid 
phase applies; we argue that it is a better strategy in developing EOS to try 
to understand the true Hamiltonian of the system, and then to use it when 
doing statistical mechanics. Almost inevitably, one must make approxima- 
tions (which we certainly have done here), but we believe we are in a better 
position to understand and improve upon them if the physical foundation of 
the EOS is as solid as we can make it. 

Finally, let us make some conservative estimates of the range of applica- 
bility of this EOS. Any limits will arise from two sources: the validity of the 
approximation that F l &h is negligible in the liquid (see the Introduction), and 
the limited ranges over which the functions §o(V), g(uj), n(e), and T m (P) 
are known. Let's consider each in turn. 

(1) We know from experiment that F^ h is negligible when T < 3 T m (again, 
see the Introduction), and judging from trends in the data we suspect Fj[ b 
will still be small up to T « 5 T m , but clearly this term must become relevant 
as the nuclear motion becomes more gaslike. Thus we shall take care with 
any data at p and T such that T approaches or exceeds 5T m (p). 

(2) At densities below approximately 6 g/cm 3 , we are confident that the 
solid is in the fee phase, and the liquid free energy is based on this phase, so 
we trust the full EOS here. At higher densities, we must be more circum- 
spect; the solid may have undergone a phase transition to bec, and the liquid 
EOS at this density may be based on the wrong solid free energy. Further, 
as we have indicated earlier, Eq. (|23|) for the melt curve has received inde- 
pendent support only up to 200 GPa, so we must be cautious with the liquid 
EOS in regions beyond this point. We decided to be brave and accept the 
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melt curve as valid up to 400 GPa; this corresponds to a liquid density of 
6.15 g/cm 3 , and since this is not far from the probable location of the solid 
fcc-bcc transition, we take it as the density limit of our EOS. (Even if we did 
not have this concern, we would be restricted to densities below 8.17 g/cm 3 , 
where electronic structure results are available.) Also, the free electron ap- 
proximation to n s (e) has been validated only for e — ep up to 1/2 Ry, or 6.8 
eV, at low compression and 1 Ry, or 13.6 eV, at high compression, but at 
higher temperatures the electronic energy and entropy are sensitive to the 
details of n s (e) to energies above these limits. We estimated the values of T 
that begin to probe the unvalidated region of n s (e) (roughly 3kT = e — ep), 
and we found that over our valid density range the T = 5T m limit always 
took precedence. Hence this limit is not relevant for us, but we mention it for 
completeness, as it may become a concern if the EOS is extended to higher 
densities. 

Figure El shows the limits p < 6.15 g/cm 3 and T < 5T m of the EOS in 
T — P space, together with the melt curve and the Hugoniot (see the next 
Subsection), while Figure [7| shows the same three features in T — p space. 
In this Figure, the melt curve becomes a two-phase region, which we will 
consider in more detail in the next Subsection. We are confident that this 
EOS is valid within these limits, but we don't know how far beyond them 
the inaccuracies begin to appear; thus we will not be shy about considering 
data not too far outside this range. 

3.2 Comparison with Hugoniot data 

If a shock wave travels at speed u s through a sample of material, accelerating 
its particles from rest to speed u p and changing its density, atomic volume, 
pressure, and internal energy per atom from p , V , P , and E to p, V, P, 
and E, then (assuming thermal equilibrium before and after the shock) these 
quantities must satisfy the Rankine-Hugoniot relations, 

p(u s - Up) = p u s 
P - Pq = PoU s u p 

E-E = ±(P + P)(V -V), (28) 

derived from considerations of mass, momentum, and energy conservation. 
(It is assumed that the wave is steady and strength effects are negligible.) 



21 




100 200 300 400 500 600 700 

P (GPa) 

Figure 6: The limits of our EOS, the melt curve, and the Hugoniot. 
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Figure 7: The limits of our EOS, the two-phase region (solid below the region, 
liquid above), and the Hugoniot. 
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Figure 8: The u s -u p Hugoniot for Al predicted by our EOS together with ex- 
perimental data from [1312011211 122 123 121 • The intersection of the Hugoniot 
with the limit of validity of the EOS (dot-dash line) is also indicated. 



24 




3 4 5 6 
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Figure 9: The P-p Hugoniot for Al predicted by our EOS together with ex- 
perimental data from [HI |HJ |2H 122 1211211 • The intersection of the Hugoniot 
with the limit of validity of the EOS (dot-dash line) is also indicated. 

By solving these equations together with the EOS, which relates P, V, 
and E, we can determine the Hugoniot, the curve of all possible end states 
of the shocked material. We used our EOS and Eqs. (J25j) to compute u s as 
a function of u p and P as a function of p along the Al Hugoniot; the results 
are shown in Figures |H1 and El along with the intersection of the Hugoniot 
with the limit of validity of the EOS. Hugoniot data from several sources 
03 1201 ED 1221 123 123 are also included. The low-pressure region of the 
Hugoniot is highlighted in Figures ITU1 and ITU and the intermediate-pressure 
region, including the intersections with the phase boundaries, is shown in 
Figures C21 and d 
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Figure 10: The u s -u p Hugoniot in the low-P region, with data from [T§1 I21[ 
221 12S] • The Up error bars on the circles |2S] appear as slightly broadened 
vertical lines. 



26 




2.8 3 3.2 3.4 3.6 3.8 4 4.2 

p (g/cm 3 ) 



Figure 11: The P-p Hugoniot in the low-P region, with data from [TT ^ I21 [ 

E2E3!. 
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Figure 12: The u s -u p Hugoniot in the intermediate- P region, including in- 
tersections with the phase boundaries, with data from \23\ 12*1] . 

Three important considerations in selecting which data to include are 
(1) the initial densities of the samples, (2) the quality of the experimental 
technique, and (3) whether the measurements were absolute or relative. All 
of the available data were taken using Al alloys with densities that differ 
from the known pure metal value of 2.70 g/cm 3 (predicted correctly by our 
EOS); some alloys are as close as 2.71 g/cm 3 while others differ much more. 
Since Hugoniots in general are quite sensitive to the initial density, we chose 
to compare only with the data for which po clustered around 2.71 g/cm 3 . 
(Thus we used only one data point from [20], which mainly concerns porous 
materials. All of the data from the other references were used.) We also 
avoided sources which gathered data using unusual shock wave geometries 
(such as PS])) an d we also chose not to use the results of indirect or relative 
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Figure 13: The P-p Hugoniot in the intermediate- P region, including inter- 
sections with the phase boundaries, with data from [23J 1211- 
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measurements, such as E3 EE] , preferring to rely on the absolute mea- 
surements that are available. Finally, we did not use the few data points 
available (primarily nuclear-driven) that lie very far beyond the limits of 
applicability of our EOS (but see below). 

The theoretical Hugoniot compares well with both the u s -u p and the P-p 
data all the way up to the predicted limit of its validity, at approximately 500 
GPa (5 Mbar). More specifically, theory agrees with experiment at P ~ 40 
GPa ( Figures ITU1 and ITT|) ; at 40-125 GPa, theory falls below the experimental 
error bars by around 1% at most (Figures [TUl through [T3J); and theory again 
lies within the experimental error bars through the liquid phase (Figures |Hl El 
H2J andHHJ). (We recall that given percentage errors in u s and u p correspond 
to roughly the same percentage errors in the P-p plane.) The presence of 
theoretical error only in the solid phase is likely due to strength effects, which 
are present in the solid but not in the liquid, and which are neglected in our 
Hugoniot calculations. Furthermore, as Figure H*3*l shows, we predict that the 
Hugoniot crosses the two-phase region between p = 4.43 g/cm 3 and p = 4.58 
g/cm 3 , corresponding to a range in P from 126 to 156 GPa; this agrees very 
well with [Hj, in which melting was found to occur between 125 and 150 
GPa. (We note that ^3] used Al 2024, an alloy whose density is sufficiently 
different from pure Al that we did not use Hugoniot data taken with that 
alloy in our Figures. We consider their melting results because, as we saw 
in the previous Subsection, their data are consistent with other experiments 
that did use pure Al.) The correction to $q from the previous Subsection, 
shown in Figure El shifts the Hugoniot at pressures below 30 GPa, bringing it 
into excellent agreement with experiment, while at pressures above 60 GPa 
or so, the effect on the Hugoniot is insignificant. 

We have also compared our results with data just beyond the EOS limits 
of validity; the points in |2H] that match our initial density (one of which is a 
reanalysis of the single point in [10]), lying at about 10 Mbar, fall noticeably 
below our Hugoniot, and their consistency with the very-high-pressure points 
of Ragan j^H E2] strongly suggests that they represent the true Hugoniot, 
which thus falls beneath our prediction at higher pressures. Possible errors 
in our EOS at such densities include, in what we estimate to be decreasing 
order of importance, (1) the shift from the fee to the bec crystal, with a 
corresponding change in $(, as discussed in the last Subsection, (2) deviations 
in the melt curve T m (P) from the Boehler-Ross form at higher pressures (the 
densities of the points in [2H] correspond to melt pressures around 620 GPa 
according to our EOS), (3) the fact that at such high T the EOS is probing 
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the high-energy region of n s (e), and (4) the neglect of the anharmonic and 
boundary contributions to the liquid EOS (T is only slightly below 5T m at 
these points according to our EOS). 

4 Conclusions 

Drawing upon theory developed in [2], we have described a framework for 
constructing EOS for elemental solids and liquids, and we have discussed 
experimental and theoretical results indicating that the framework remains 
highly accurate at low pressures when certain small effects (anharmonicity, 
boundaries, electron-phonon interactions) are neglected. After displaying 
the resulting formulas for the Helmholtz free energy, we considered the in- 
formation one needs to evaluate them, and we discussed the combination of 
experiment and theoretical work that could be used to get this information. 
Finally, as an illustration we constructed an EOS for Al, established its range 
of validity based on the inputs to the EOS, and compared it with Hugoniot 
data to 5 Mbar; our EOS matched the data to the accuracy we expected 
based on the low-pressure results. 

We consider the primary advantage of this method to lie in the fact that 
it incorporates into the decomposition of the Hamiltonian a great deal of 
accumulated knowledge of condensed matter physics both for the solid and 
liquid phases (for example, the fact that the electronic ground state energy 
is the most appropriate potential for the nuclear motion). If we have indeed 
captured the correct physics (and we expect no new physics to enter until 
the relativistic domain), then the EOS should have the right functional form, 
which means that if it is shown to agree with available data, then we have 
reason to believe that it will be equally accurate in regions where no data 
are available; and making predictions where we have no data is the point of 
having an EOS to begin with. Furthermore, the better the foundation we can 
build, the better our position for intelligently investigating and controlling 
our approximations. 

This discussion bears on the second goal of this paper, which was to 
learn whether the approximation of neglecting anharmonic, boundary, and 
electron-phonon effects remains useful at higher pressures. We already knew, 
as discussed in the Introduction, that at low P the anharmonic and electron- 
phonon terms are small, and we found this from direct calculations; we also 
knew that for several elements, over a range of T, at low P the approxi- 
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mations in question yielded thermal energies and entropies that disagreed 
with experiments by 5% at most. Our work here has shown that for one 
material at much higher P the approximations yield results that match data 
along a single curve, the Hugoniot, to comparable accuracy. Based on our 
arguments above, that the EOS incorporates the correct physics and is thus 
of the correct functional form, we claim to have shown that this Al EOS is 
trustworthy throughout its range of validity, for all T and P. 

The main disadvantage of this method is that it relies on many inputs 
($o, g{u), and n{e) for each phase) which may be available only over limited 
ranges, and each of these limits also restricts the range of validity of the EOS. 
Our Al example amply illustrates this problem; with a compression range 
from a little under one to just over two, and a temperature range that reaches 
only to slightly under 4 eV, this EOS is inadequate for many applications 
at the national laboratories. We argue, however, that this problem does 
not indicate a deficiency in the approach; it only underscores the need for 
many more DFT calculations of these quantities for more materials with ever 
greater accuracy over ever larger ranges. 

In the meantime, though, we would like to be able to say something about 
elemental solid and liquid EOS at higher compressions. We do know that as 
density increases, the electrons come to dominate the free energy, and it is 
also known that TFD correctly describes the electrons in the limit of infinite 
density. This suggests the following possibility: Construct an EOS using the 
present techniques to compressions as high as the available experimental or 
DFT results allow, and then interpolate between these results and the pre- 
dictions of TFD for higher compressions. This raises an important question: 
At what pressures does TFD begin to become accurate? Conventional wis- 
dom, usually traced back to Feynman et al. [33] , has held that TFD becomes 
reliable starting at P « 10 Mbar, but other work 34 a suggests that TFD 
(or TF in their case, but TF and TFD converge at high pressures) deviates 
noticeably from electronic structure results until 100 Mbar at least. This 
suggests that the pressure threshold at which TFD is trustworthy has not 
yet been adequately established; it would be of great interest to settle this 
question more definitively. 
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